
#include "float.h"
#define blockSize 256





void check_cache(int iHigh, int iLow, bool& iHighCompute, bool& iLowCompute, int& iHighCacheIndex, int& iLowCacheIndex, int* Cache, int cacheSize){


	iHighCacheIndex = iHigh % cacheSize;
	iLowCacheIndex = iLow % cacheSize;

	//iHigh is the cache
	if(Cache[iHighCacheIndex] == iHigh) iHighCompute = false;
	else iHighCompute = true;

	if(Cache[iLowCacheIndex] == iLow) iLowCompute = false;
	else iLowCompute = true;

	//When updating the cache, we need to do the iHigh part first. Because when iHighCacheIndex==iLowCacheIndex, iLow will overwrite iHigh position.	
	if(iHighCompute) Cache[iHighCacheIndex] = iHigh;
	if(iLowCompute) Cache[iLowCacheIndex] = iLow;

}




static char *line = NULL;
static int max_line_len;


static char* readline(FILE *input)
{
	int len;

	if(fgets(line,max_line_len,input) == NULL)
		return NULL;

	while(strrchr(line,'\n') == NULL)
	{
		max_line_len *= 2;
		line = (char *) realloc(line,max_line_len);
		len = (int) strlen(line);
		if(fgets(line+len,max_line_len-len,input) == NULL)
			break;
	}
	return line;
}


int readSvm(const char* filename, float** p_data, float** p_labels, int* p_nObs, int* p_nDim) {

	FILE* inputFile = fopen(filename, "r");
	if (inputFile == 0) {
		printf("File not found\n");
		return 0;
	}

	int npoints, dimension,inst_max_index, index;
	float value;
	char *endptr;
	char *idx, *val, *label, *p;

	npoints = 0;
	dimension = 0;
	inst_max_index = 0;
	max_line_len = 1024;
	line = (char*) malloc(max_line_len * sizeof(char));

	//determine the number of training samples and the number of features
	while(readline(inputFile)!=NULL)
	{
		p = strtok(line," \t");
		while(1)
		{

			idx = strtok(NULL,":");
			p = strtok(NULL," \t");

			if(p == NULL || *p == '\n') // check '\n' as ' ' may be after the last feature
				break;

			inst_max_index = (int) strtol(idx,&endptr,10);

			if(inst_max_index > dimension)
				dimension = inst_max_index;


			
		}
		++npoints;
	}

	rewind(inputFile);
	*(p_nObs) = npoints;
	*(p_nDim) = dimension;


	float* data = (float*)malloc(sizeof(float)*npoints*dimension);
	float* labels = (float*)malloc(sizeof(float)*npoints);
	*(p_data) = data;
	*(p_labels) = labels;


	for(int i = 0; i < npoints; i++) {

		readline(inputFile);
		label = strtok(line," \t\n");
		labels[i] = atof(label);


		while(1)
		{
			idx = strtok(NULL,":");
			val = strtok(NULL," \t");

			if(val == NULL)
				break;

			index = (int) strtol(idx,&endptr,10);
			value = (float) strtod(val,&endptr);

			data[i*dimension + index -1] =value;
		}
	}

	fclose(inputFile);
	return 1;
}

void printModel(const char* outputFileName, float* alpha, float* labels, float* data, int nPoints, int nDimension, float epsilon, float rho) { 

	printf("Output File: %s\n", outputFileName);
	FILE* outputFile = fopen(outputFileName, "w");
	if (outputFile == NULL) {
		printf("Can't write %s\n", outputFile);
		exit(1);
	}

	int nSV = 0;
	int pSV = 0;
	for(int i = 0; i < nPoints; i++) {
		if (alpha[i] > epsilon) {
			if (labels[i] > 0) {
				pSV++;
			} else {
				nSV++;
			}
		}
	}


	fprintf(outputFile,"svm_type c_svc\n");
	fprintf(outputFile,"kernel_type linear\n");
	fprintf(outputFile, "nr_class 2\n");
	fprintf(outputFile, "total_sv %d\n",pSV+nSV);
	fprintf(outputFile, "rho %.10f\n", rho);
	fprintf(outputFile, "label 1 -1\n");
	fprintf(outputFile, "nr_sv %d %d\n", pSV, nSV);
	fprintf(outputFile, "SV\n");

	for (int i = 0; i < nPoints; i++) {
		if (alpha[i] > epsilon) {
			fprintf(outputFile, "%f ", labels[i]*alpha[i]);
			for (int j = 0; j < nDimension; j++) {
				fprintf(outputFile, "%d:%f ", j+1, data[j*nPoints + i]);
			}
			fprintf(outputFile, "\n");
		}
	}
	fclose(outputFile);
}




///////////////////////////// LINEAR KERNEL FUNCTIONS ////////////////////////////////
//compute phi(A'B)
static __device__ __host__ float kernel(float *pA, float *pB, int nDim) {
	float acc = 0.0f;
	for(int i=0; i<nDim; i++)
		acc += (*(pA+i)) * (*(pB+i));
	return acc;
}
//compute phi(a,b) and phi(a,c), where b and c are already in shared memory
//and a is in global memory
static __device__ __host__ 
	void dualKernel(float *pA, float *pB, float *pC, int nDim, float& phi1, float& phi2) {
		float acc1 = 0.0f; float acc2 = 0.0f;
		for(int i=0; i<nDim; i++){
			float tmpA = *(pA+i);
			acc1 += tmpA * (*(pB+i));
			acc2 += tmpA * (*(pC+i));
		} 
		phi1 = acc1; phi2 = acc2;
}

////////////////////////////// REDUCTION KERNELS //////////////////////////////////////
__device__ void argMin(int AIdx, float AVal, int BIdx, float BVal, int* idx, float* val) {
	if (BVal < AVal) { *idx = BIdx; *val = BVal; }
	else { *idx = AIdx; *val = AVal; }
}
__device__ void argMax(int AIdx, float AVal, int BIdx, float BVal, int* idx, float* val) {
	if (BVal > AVal) { *idx = BIdx; *val = BVal; }
	else { *idx = AIdx; *val = AVal; }
}

	

__device__ void ReduceMin(float* values, int* indices){
	int tid = threadIdx.x;
	if(blockSize >=512){if(tid<256){argMin(indices[tid], values[tid], indices[tid+256], values[tid+256], indices+tid, values+tid);} __syncthreads();}
if(blockSize >=256){if(tid<128){argMin(indices[tid], values[tid], indices[tid+128], values[tid+128], indices+tid, values+tid);} __syncthreads();}
if(blockSize >=128){if(tid< 64){argMin(indices[tid], values[tid], indices[tid+ 64], values[tid+ 64], indices+tid, values+tid);} __syncthreads();}
if(tid<32){
	if(blockSize>=64)  argMin(indices[tid], values[tid], indices[tid+32], values[tid+32], indices+tid, values+tid);
	if(blockSize>=32)  argMin(indices[tid], values[tid], indices[tid+16], values[tid+16], indices+tid, values+tid);
	if(blockSize>=16)  argMin(indices[tid], values[tid], indices[tid+ 8], values[tid+ 8], indices+tid, values+tid);
	if(blockSize>= 8)  argMin(indices[tid], values[tid], indices[tid+ 4], values[tid+ 4], indices+tid, values+tid);
	if(blockSize>= 4)  argMin(indices[tid], values[tid], indices[tid+ 2], values[tid+ 2], indices+tid, values+tid);
	if(blockSize>= 2)  argMin(indices[tid], values[tid], indices[tid+ 1], values[tid+ 1], indices+tid, values+tid);
}
}



__device__ void ReduceMax(float* values, int* indices){
	int tid = threadIdx.x;
	if(blockSize >=512){if(tid<256){argMax(indices[tid], values[tid], indices[tid+256], values[tid+256], indices+tid, values+tid);} __syncthreads();}
if(blockSize >=256){if(tid<128){argMax(indices[tid], values[tid], indices[tid+128], values[tid+128], indices+tid, values+tid);} __syncthreads();}
if(blockSize >=128){if(tid< 64){argMax(indices[tid], values[tid], indices[tid+ 64], values[tid+ 64], indices+tid, values+tid);} __syncthreads();}
if(tid<32) {
	if(blockSize>=64)  argMax(indices[tid], values[tid], indices[tid+32], values[tid+32], indices+tid, values+tid);
	if(blockSize>=32)  argMax(indices[tid], values[tid], indices[tid+16], values[tid+16], indices+tid, values+tid);
	if(blockSize>=16)  argMax(indices[tid], values[tid], indices[tid+ 8], values[tid+ 8], indices+tid, values+tid);
	if(blockSize>= 8)  argMax(indices[tid], values[tid], indices[tid+ 4], values[tid+ 4], indices+tid, values+tid);
	if(blockSize>= 4)  argMax(indices[tid], values[tid], indices[tid+ 2], values[tid+ 2], indices+tid, values+tid);
	if(blockSize>= 2)  argMax(indices[tid], values[tid], indices[tid+ 1], values[tid+ 1], indices+tid, values+tid);
}
}


/////////////// SEQUENTIAL MINIMAL OPTIMIZATION KERNELS ///////////////////////////////
__global__ void smo_init(float* dData, int dDataPitchSize, int nObs, int nDim, float* dKernelDiag, float* dAlpha, float* dF, float* dLabels) { 
	int idx = blockIdx.x*blockDim.x + threadIdx.x;

	//compute diagonal elements of kernel matrix (each thread does a different one)
	//compute its own F (in this case, initialize to the value of the label)
	//set alphas to 0
	if (idx < nObs) {
		dKernelDiag[idx] = kernel(dData+idx*dDataPitchSize,dData+idx*dDataPitchSize, nDim);
		dF[idx] = -dLabels[idx];
		//dF[idx] = -1;
		dAlpha[idx] = 0;
	}
}

__global__ void smo_first_iter(float* dResult, float* dKernelDiag, float* dData, int dDataPitchSize, float* dAlpha, float cost, int nDim, int iLow, int iHigh) { 

	float* pA = dData + iHigh * dDataPitchSize;
	float* pB = dData + iLow * dDataPitchSize;
	float eta = dKernelDiag[iHigh] + dKernelDiag[iLow] - 2*kernel(pA,pB,nDim);
	float alpha2New = 2/eta;
	if (alpha2New > cost) alpha2New = cost;
	if (alpha2New<0) alpha2New =0;
	dAlpha[iLow] = alpha2New;
	dAlpha[iHigh] = alpha2New;

	//return results    
	*(dResult + 0) = 0.0;
	*(dResult + 1) = 0.0;
	*(dResult + 2) = 1.0;
	*(dResult + 3) = -1.0;
	*(dResult + 6) = alpha2New;
	*(dResult + 7) = alpha2New;
}


__global__ void smo_local(float* dData, int dDataPitchSize,  float* dLabels, int nObs, int nDim, float epsilon, float cost_epsilon, float* dAlpha, float* dF, float alpha1Diff, float alpha2Diff, int iLow, int iHigh, float* dCache, int dCachePitchSize, int iLowCacheIndex, int iHighCacheIndex, int* devLocalIndicesRL, int* devLocalIndicesRH, float* devLocalFsRL, float* devLocalFsRH, bool iHighCompute, bool iLowCompute){

	//xIHighLow size = 2*nDim
	extern __shared__ float xIHighLow[];
	__shared__ int tempLocalIndices[blockSize];
	__shared__ float tempLocalFs[blockSize];
    int globalIndex = blockDim.x*blockIdx.x + threadIdx.x;

	if (iHighCompute) {
		//Load xIHigh into shared memory
		if(threadIdx.x==0)
		{
			for(int i =0; i< nDim; i++)
				xIHighLow[i] = dData[iHigh*dDataPitchSize + i];
		}
	}

	if (iLowCompute) {
		//Load xILow into shared memory
		if(threadIdx.x==0)
		{
			for(int i=0; i<nDim; i++)
				xIHighLow[nDim+i]=dData[iLow*dDataPitchSize + i];
        }
	}

	__syncthreads();

	float alpha;
	float f;
	float label;
	int reduceFlag; //ihigh =0, ilow=1, both=2, none=-1

	if (globalIndex < nObs) {
		alpha = dAlpha[globalIndex];
		f = dF[globalIndex];
		label = dLabels[globalIndex];

		if (alpha > epsilon) {
			if (alpha < cost_epsilon) {
				reduceFlag = 2; //Unbound support vector (I0)
			} else {
				if (label > 0) {
					reduceFlag = 1; //Bound positive support vector (I3)
				} else {
					reduceFlag = 0; //Bound negative support vector (I2)
				}
			}
		} else {
			if (label > 0) {
				reduceFlag = 0; //Positive nonsupport vector (I1)
			} else {
				reduceFlag = 1; //Negative nonsupport vector (I4)
			}
		}
	} else {
		reduceFlag = -1;
	}
	
		float highKernel = 0;
		float lowKernel = 0;
		if (reduceFlag >= 0) {
			if (!iHighCompute) {
				highKernel = dCache[(dCachePitchSize * iHighCacheIndex) + globalIndex];
			}
			if (!iLowCompute) {
				lowKernel = dCache[(dCachePitchSize * iLowCacheIndex) + globalIndex];
			}
			if (iHighCompute && iLowCompute) {
				dualKernel(dData + globalIndex*dDataPitchSize, xIHighLow, xIHighLow + nDim ,nDim, highKernel, lowKernel);
			} else if (iHighCompute) {
				highKernel = kernel(dData + globalIndex * dDataPitchSize, xIHighLow, nDim);
			} else if (iLowCompute) {
				lowKernel = kernel(dData + globalIndex * dDataPitchSize,  xIHighLow + nDim, nDim);
			}

			//highKernel = phi(x_i, x_high);
			//lowKernel = phi(x_i, x_low);
			//Update the f_i
			//f = f + alpha1Diff * highKernel * dLabels[globalIndex];
			//f = f + alpha2Diff * lowKernel * dLabels[globalIndex];

			f = f + alpha1Diff * highKernel;
			f = f + alpha2Diff * lowKernel;

			//update the cache
			if (iHighCompute) dCache[(dCachePitchSize * iHighCacheIndex) + globalIndex] = highKernel;
			if (iLowCompute)  dCache[(dCachePitchSize * iLowCacheIndex) + globalIndex] = lowKernel;
				
			dF[globalIndex] = f;

		}
		__syncthreads();

		//Ilow
		if ((reduceFlag==1) || (reduceFlag==2)) {
			//tempLocalFs[threadIdx.x] = -f * dLabels[globalIndex];
			tempLocalFs[threadIdx.x] = f;
			tempLocalIndices[threadIdx.x] = globalIndex;
		} else {
			tempLocalFs[threadIdx.x] = -FLT_MAX; //Ignore me
		}

		__syncthreads();
		ReduceMax(tempLocalFs, tempLocalIndices);
		//argmaxReduce(tempLocalFs, tempLocalIndices);

		__syncthreads();


		if (threadIdx.x == 0) {
			devLocalIndicesRL[blockIdx.x] = tempLocalIndices[0];
			devLocalFsRL[blockIdx.x] = tempLocalFs[0];
		}

		__syncthreads();

		//Ihigh
		if ((reduceFlag ==0)|| (reduceFlag ==2)) {
			//tempLocalFs[threadIdx.x] = -f * dLabels[globalIndex];
			tempLocalFs[threadIdx.x] = f;
			tempLocalIndices[threadIdx.x] = globalIndex;
		} else {
			tempLocalFs[threadIdx.x] = FLT_MAX; //Ignore me
		}
		__syncthreads();
		ReduceMin(tempLocalFs, tempLocalIndices);
		//argminReduce(tempLocalFs, tempLocalIndices);
		if (threadIdx.x == 0) {
			devLocalIndicesRH[blockIdx.x] = tempLocalIndices[0];
			devLocalFsRH[blockIdx.x] = tempLocalFs[0];
		}
}


__global__ void smo_global(float* dData, int dDataPitchSize, float* dLabels, float* dKernelDiag, float* dAlpha, float* dResult, float cost, int nDim, int* devLocalIndicesRL, int* devLocalIndicesRH, float* devLocalFsRL, float* devLocalFsRH, int inputSize, float epsilon){


	__shared__ int tempIndices[blockSize];
	__shared__ float tempFs[blockSize];

	int iHigh;
	int iLow;
	float bHigh;
	float bLow;

	//Load elements
	if (threadIdx.x < inputSize) {
		tempIndices[threadIdx.x] = devLocalIndicesRH[threadIdx.x];
		tempFs[threadIdx.x] = devLocalFsRH[threadIdx.x];
	} else {
		tempFs[threadIdx.x] = FLT_MAX;
	}


	if (inputSize > blockSize) {
		for (int i = threadIdx.x + blockSize; i < inputSize; i += blockDim.x) {
			argMin(tempIndices[threadIdx.x], tempFs[threadIdx.x], devLocalIndicesRH[i], devLocalFsRH[i], tempIndices + threadIdx.x, tempFs + threadIdx.x);
		}
	}
	__syncthreads();
	//argminReduce(tempFs, tempIndices);
	ReduceMin(tempFs, tempIndices);

	__syncthreads();
	if(threadIdx.x==0){
		iHigh = tempIndices[0];
		bHigh = tempFs[0];
	}

	__syncthreads();

	//Load elements
	if (threadIdx.x < inputSize) {
		tempIndices[threadIdx.x] = devLocalIndicesRL[threadIdx.x];
		tempFs[threadIdx.x] = devLocalFsRL[threadIdx.x];
	} else {
		tempFs[threadIdx.x] = -FLT_MAX;
	}


	//blockDim.x==number per block
	if (inputSize > blockSize) {
		for (int i = threadIdx.x + blockSize; i < inputSize; i += blockDim.x) {
			argMax(tempIndices[threadIdx.x], tempFs[threadIdx.x], devLocalIndicesRL[i], devLocalFsRL[i], tempIndices + threadIdx.x, tempFs + threadIdx.x);
		}
	}
	__syncthreads();

	ReduceMax(tempFs, tempIndices);
	//argmaxReduce(tempFs, tempIndices);

	__syncthreads();

	if (threadIdx.x == 0) {

		iLow = tempIndices[0];
		bLow = tempFs[0];

		float eta = dKernelDiag[iHigh] + dKernelDiag[iLow];

		float kernelEval = kernel(dData+iHigh*dDataPitchSize, dData + iLow *dDataPitchSize, nDim);

		eta = eta - 2*kernelEval;
		// clipping method from libsvm
		//int ylow = dLabels[iLow];
		//int yhigh = dLabels[iHigh];
/*
		eta = eta - 2*ylow*yhigh*kernelEval;
		if(eta<=0) eta = epsilon;

		float alpha1Old = dAlpha[iHigh];
		float alpha2Old = dAlpha[iLow];

		float b = -ylow * bLow + yhigh*bHigh;

		float alpha2New = alpha2Old + ylow * b/eta;
		float alpha1New = alpha1Old - yhigh *b/eta;

		float sum = ylow * alpha2Old + yhigh * alpha1Old;

		if(alpha2New > cost) alpha2New = cost;
		if(alpha2New < 0) alpha2New =0;

		alpha1New = yhigh*(sum - ylow*alpha2New);
		
		if(alpha1New > cost) alpha1New = cost;
		if(alpha1New < 0) alpha1New =0;

		alpha2New = ylow*(sum - yhigh*alpha1New);

		float alpha1Old = dAlpha[iHigh];
		float alpha2Old = dAlpha[iLow];

		int ylow = dLabels[iLow];
		int yhigh = dLabels[iHigh];

		float alpha2New = alpha2Old + ylow*(bHigh - bLow) / eta;

		if(alpha2New<0) alpha2New =0;
		if(alpha2New>cost) alpha2New = cost;

		float alpha1New = alpha1Old + ylow * yhigh * (alpha2Old - alpha2New);
		if(alpha1New<0) alpha1New =0;
		if(alpha1New>cost) alpha1New = cost;

*/
		float alpha1Old = dAlpha[iHigh];
		float alpha2Old = dAlpha[iLow];
		float alphaDiff = alpha2Old - alpha1Old;
		float lowLabel = dLabels[iLow];
		float sign = dLabels[iHigh] * lowLabel;
		float alpha2UpperBound;
		float alpha2LowerBound;
		if (sign < 0) {
			if (alphaDiff < 0) {
				alpha2LowerBound = 0;
				alpha2UpperBound = cost + alphaDiff;
			} else {
				alpha2LowerBound = alphaDiff;
				alpha2UpperBound = cost;
			}
		} else {
			float alphaSum = alpha2Old + alpha1Old;
			if (alphaSum < cost) {
				alpha2UpperBound = alphaSum;
				alpha2LowerBound = 0;
			} else {
				alpha2LowerBound = alphaSum - cost;
				alpha2UpperBound = cost;
			}
		}
		float alpha2New;
		if (eta > 0) {
			alpha2New = alpha2Old + lowLabel*(bHigh - bLow)/eta;
			if (alpha2New < alpha2LowerBound) {
				alpha2New = alpha2LowerBound;
			} else if (alpha2New > alpha2UpperBound) {
				alpha2New = alpha2UpperBound;
			}
		} else {
			float slope = lowLabel * (bHigh - bLow);
			float delta = slope * (alpha2UpperBound - alpha2LowerBound);
			if (delta > 0) {
				if (slope > 0) {
					alpha2New = alpha2UpperBound;
				} else {
					alpha2New = alpha2LowerBound;
				}
			} else {
				alpha2New = alpha2Old;
			}
		}
		float alpha2Diff = alpha2New - alpha2Old;
		float alpha1Diff = -sign*alpha2Diff;
		float alpha1New = alpha1Old + alpha1Diff;



		*(dResult + 0) = alpha2Old;
		*(dResult + 1) = alpha1Old;
		*(dResult + 2) = bLow;
		*(dResult + 3) = bHigh;
		dAlpha[iLow] = alpha2New;
		dAlpha[iHigh] = alpha1New;
		*(dResult + 4) = alpha2New;
		*(dResult + 5) = alpha1New;
		*(dResult + 6) = float(iLow);
		*(dResult + 7) = float(iHigh);
	}

}

